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DFT 


My  efforts  in  the  Density  Functional  Theory  (DFT)  project  of  Dr. 
Bartlett's  research  group  have  been  limited  to  the  implementation  of 
analytical  energy  gradients  with  respect  to  nuclear  displacement  in  ACES2. 
The  most  immediate  benefit  of  this  work  is  analytical  geometry  optimizations 
and  finite-difference  vibrational  frequencies  with  analytical  gradients. 

In  the  most  general  case,  the  Kohn-Sham  (KS)  energy  equation  looks 
almost  identical  to  the  well-known  Hartree-Fock  (HF)  energy  expression. 

This  is  not  surprising  since  the  HF  non-local  exchange  functional,  -K,  is  one 
of  many  possible  KS  exchange  operators.  This  allows  for  a  relatively  simple 
"port"  of  a  pre-existing  HF  Self-Consistent  Field  (SCF)  program  to  a  more 
robust  KS  SCF  program. 

DFT  was  added  to  ACES  many  years  ago  in  the  form  of  tack-on 
functionals,  also  known  as  "Hartree-Fock  Density  Functional  Theory,"  in 
which  the  HF  density  was  used  to  generate  energies  with  well-known  energy 
functionals  like  BLYP,  B3LYP,  and  LDA-VWN.  Shortly  thereafter,  KS  DFT  was 
added  to  ACES  as  a  complete  re-write  of  the  existing  HF  SCF  program.  The 
robust  numerical  integrator  accompanying  this  program  forms  the 
cornerstone  of  the  KS  DFT  analytical  gradients. 


KS  DFT  analytical  gradients 

The  HF  energy  gradient  is  given  by  equation  1.  The  Greek  subscripts 
denote  atomic  orbital  basis  functions  and  the  capital  superscripts  denote 
differentiation  with  respect  to  that  nuclear  center.  The  P,  h,  W,  and  S 
matrices  are  the  density,  one-particle  operator,  energy-weighted  density, 
and  AO  overlap  matrices  respectively. 
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The  KS  energy  gradient  is  given  by  equation  2.  The  most  notable 
change  is  the  separation  of  the  double-bar  integral  from  (l-P(12))/ri2  to 
l/ri2  +  Vxc.  Choosing  Vxc  to  be  the  HF  non-local  exchange  operator,  -K, 
readily  reduces  eqn.  2  to  eqn.  1. 

i;E = -  5XX 

V  ™  J  M’v  (2) 

+  E^v((^Kc|v)  +  (^K|va)) 

J U,V 

Adding  analytical  KS  energy  gradients  to  ACES  was  a  two-step 
process:  1)  parameterize  the  -K  contribution  in  the  existing  gradient  code 
and  2)  add  the  exchange-correlation  contribution  to  the  gradient  as  a  stand¬ 
alone  program.  Both  of  these  were  completed  and  the  qualities  of  the  results 
are  mixed. 

We  do  not  incorporate  the  derivative  of  the  integration  grid  with 
respect  to  each  nuclear  center  into  the  energy  gradient.  This  immediately 
introduces  an  error,  the  size  of  which  depends  upon  the  granularity  of  the 
grid.  Using  the  default  grid  for  both  calculations,  the  difference  between  the 
analytically-  and  numerically-optimized  geometries*  of  H2  is  approximately 
lxlO'4  Angstroms.  Although  this  seems  large,  the  actual  energy  difference  is 
lxlO'8  a.u. 

I  claim  these  results  are  of  mixed  quality  because  for  larger  systems, 
those  with  more  geometric  parameters  to  optimize  and  those  of  lower 
molecular  symmetry,  the  agreement  between  numerically-  and  analytically- 
optimized  geometries  goes  down;  however,  the  energy  differences  agree 
well-within  chemical  accuracy.  The  problem  with  this  effect  is  that  the  user 
of  the  program  will  immediately  begin  to  suspect  the  correctness  of  the 
results,  especially  when  compared  to  other  DFT  programs,  when  the  energy 


*  The  optimizations  and  final  energies  of  H2  were  done  in  the  DZP  basis  using  the 
B3LYP  potential  and  functional. 
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is  actually  minimized  but  the  nuclear  coordinates  do  not  appear  to  coincide 
with  other  published  data. 

The  next  two  phases  of  this  project  would  be  adding  analytical  energy 
hessians  and  integration  grid  differentiation. 

Ethylene 

One  of  the  fundamental  strengths  of  Kohn-Sham  DFT  lies  in  the  fact 
that  the  orbital  eigenvalue  of  the  HOMO  is  theoretically  the  exact  ionization 
potential  (IP)  due  to  the  correct  long-range  behavior  of  the  density.  We 
briefly  examined  the  entire  spectrum  of  eigenvalues  returned  by  KS  DFT  in 
an  attempt  to  determine  whether  the  other  orbital  energies  meant  as  much. 

The  prototypical  molecule  C2H4  was  chosen  as  a  test  system  for  its 
high  symmetry  and  low  number  of  electrons.  Of  its  eight  IPs,  the  first  five 
are  all  of  different  symmetry  giving  us  a  convenient  way  of  cataloging  the 
orbital  energies.  Not  only  were  the  low-lying  orbital  energies  off  by  an 
average  of  0.75  eV  from  experiment  but  also  the  first  principle  IP,  which 
should  be  exact,  was  off  by  0.7  eV.  This  lead  us  to  question  whether  our 
current  set  of  highly  correlated  methods  could  more  accurately  predict  the 
IPs. 

Our  most  advanced  treatment  of  IPs  is  IP-EOM-CCSD.  Using  the 
Equation-of-Motion  on  a  CCSD  single-reference  state,  one  calculation  is  able 
to  return  the  entire  spectrum  of  excitation  energies  by  a  matrix 
diagonalization.  Results  were  obtained  for  both  HF  and  KS  orbitals  serving  as 
the  unperturbed  ground  states  for  the  CCSD  calculation.  The  differences 
were  negligible  if  any,  e.g.,  the  largest  variation  was  for  the  principle  IP  and 
the  energy  difference  between  the  two  references  was  only  0.005  eV. 

We  considered  both  the  adiabatic  and  vertical  IPs  in  the  study+. 
Experimental  results  show  ethylene  cations  only  adiabatically  relax  for  the 


f  Currently  in  press. 
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first  two  ionized  states.  Zero-point  corrected  adiabatic  IPs  for  these  two 
states  obtained  with  IP-EOM-CCSD  differed  less  than  0.05  eV  from 
experiment. 

The  vertical  IPs,  defined  to  be  the  energy  difference  between  the 
neutral  and  cation  both  at  the  optimum  neutral  geometry,  differed  from  0.07 
eV  to  0.28  eV.  Our  target  error  was  0.1  eV.  The  only  way  our  best  method 
can  be  enhanced  is  by  the  inclusion  of  a  Franck-Condon  analysis,  which 
weighs  the  most  probable  vibrational  transitions  and  returns  a  distribution  of 
IPs  for  each  observed  band.  In  any  event,  the  quality  of  the  KS  orbital 
eigenvalues  is  quite  low  compared  to  ab  initio  treatments  of  correlation. 

ACES2 

My  involvement  with  the  underlying  source  code  of  ACES2  has 
contributed  to  other  efforts  as  well.  The  limited  integral  program  currently 
employed  was  replaced  by  a  more  robust  package  that  will  also  speed  the 
implementation  of  direct  integrals  into  our  SCF  treatment.  Also,  the  ACES2 
program  system  was  officially  ported  to  the  set  of  GNU  compilers  that  are 
freely  available  for  such  popular  Unices  and  Unix-like  operating  systems  as 
FreeBSD,  Solaris,  and  Linux. 

MOLCAS  integrals 

Recently,  in  collaboration  with  Drs.  Roland  Lindh  and  Ajith  Perera  and 
Mr.  Ken  Wilson,  the  aging  integral  program,  VMOL,  was  replaced  by  the 
more  robust  and  faster  integral  program  of  the  MOLCAS  suite  named 
Seward.  The  two  main  advantages  of  this  replacement  are  1)  near-unlimited 
numbers  of  AO  basis  functions  are  now  allowed  in  the  SCF  treatment  and  2) 
Seward  was  programmed  with  the  most  modern  techniques  and  is  currently 
supported  by  the  MOLCAS  developers  (unlike  VMOL  which  is  over  20  years 
old  and  resides  in  the  public  domain). 


4 


The  255  basis  function  limit  was  imposed  by  VMOL's  indexing  scheme 
wherein  all  four  indices  of  a  two-electron  integral  were  packed  into  1  32-bit 
integer.  On  systems  with  native  64-bit  integer  support  and  capable 
compilers,  this  scheme  was  trivially  expanded  to,  in  principle,  allow  up  to 
65535  AO  basis  functions.  The  new  method  for  storing  integrals  bypasses 
indices  altogether  and  implicitly  sorts  the  four-index  integrals  in  petite  lists 
so  there  is  no  limit  on  the  size  of  the  basis  set  regardless  of  architecture. 

Seward  also  has  the  capability  to  be  used  as  a  library  of  routines  for 
calculating  two-electron  integrals  directly.  This  is  opposite  the  way  in  which 
we  are  currently  using  Seward,  which  is  as  a  stand-alone  integral  file 
generator.  The  next  logical  extension  of  our  work  will  be  to  add  direct 
integral  capabilities  to  the  SCF  program  using  this  library. 

I  must  also  mention  the  fact  that  our  Coupled-Cluster  code  does  not 
use  explicit  indices  for  storing  the  T  amplitudes  but  rather  implicit  sorting. 
Therefore,  with  the  proper  integral  transformations,  CC  calculations  in  ACES 
will  also  be  able  to  take  advantage  of  basis  sets  with  more  than  255 
functions. 

Complete  GNU  port 

As  consumer-level  PCs  become  more  powerful,  networks  of  commodity 
off-the-shelf  computers  are  taking  their  places  among  the  ranks  of  central 
servers  and  supercomputers.  "Distributed  computing"  is  the  fad  and  Linux 
and  Beowulf  are  the  buzzwords.  Any  computer  or  network  of  computers 
running  Linux  will  have  the  GNU  utilities  and  compilers.  They  are  free  under 
the  GNU  General  Public  License  and  have  been  ported  to  many,  many  other 
operating  systems. 

To  date,  however,  the  ACES  system  was  never  fully  ported,  i.e.,  built 
and  tested,  to  this  set  of  tools.  Lack  of  familiarity  with  the  tools'  nuances  and 
lack  of  computing  power  from  the  machines  running  the  tools  kept  the 
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previous  ACES  developers  from  investing  any  time  in  it.  Given  the 
circumstances  of  today,  it  only  seemed  natural  to  finish  the  port. 

Now,  the  ACES2  program  is  officially  ported  to  the  GNU  utilities 
running  on  any  operating  system.  The  binaries  were  checked  for  numerical 
accuracy  and  all  of  them  were  successful.  Our  group  has  a  PC  running 
FreeBSD  with  384MB  of  memory  and  20GB  of  scratch  disk  space.  Crude 
initial  benchmarks  show  this  machine  is  more  than  capable  of  handling 
modest  test  systems.  The  same  systems  executed  on  our  Sun  Microsystems 
Enterprise  5000  server  had  only  a  25%  speed-up.  The  price  difference, 
however,  is  over  140%. 

The  advanced  capabilities  of  our  program  can  now  extend  to  research 
groups  with  less  funding,  and  our  program  can  be  installed  on  the  most 
basic  off-the-shelf  PC  server. 
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